# Estimate the models with unconditional effects of elections and the conditional effects of populist involvement
# and polarization in the main text.




#############
# load data #
#############

library(brms)
library(ggplot2)
library(ggpubr)
library(stargazer)
library(ggrepel)
library(dplyr)

load("PSRM Replication Files/PresidentialElectionsPosts.RData")

colnames(president_df)

all(c("populist_present", "polarization") %in% colnames(president_df)) # TRUE

length(unique(president_df$country)) # 26
length(unique(president_df$electionname)) # 29
length(unique(president_df$party_country)) # 52

# focus on 15 days before and after the election
day_range <- 15

sub <- president_df[which(president_df$daysinceelection >= -1*day_range
                          & president_df$daysinceelection <= day_range),]

# winners and losers
winners <- sub[which(sub$win == 1),]
losers <- sub[which(sub$win != 1),]

length(unique(winners$electionname[which(!is.na(winners$polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$economic_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$social_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_present))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_dummy))])) # 27

length(unique(losers$electionname[which(!is.na(losers$polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$economic_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$social_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_present))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_dummy))])) # 27




##################
# winners + love #
##################

# table M.1 model 1
winner_love1 <- brm(loveprop ~ post
                    + populist_present + polarization
                    + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                    + Xl + Xr 
                    + (1 + Xl + Xr | party_election),
                    data=winners,
                    warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(winner_love1)[[14]][,5])
#plot(winner_love1)
#summary(winner_love1)

# table M.2 model 1
winner_love2 <- brm(loveprop ~ post*populist_present + post*polarization
                    + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                    + Xl + Xr 
                    + (1 + Xl + Xr | party_election),
                    data=winners,
                    warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(winner_love2)[[14]][,5])
#plot(winner_love2)
#summary(winner_love2)




###################
# winners + angry #
###################

# table M.1 model 2
winner_angry1 <- brm(angryprop ~ post
                     + populist_present + polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=winners,
                     warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(winner_angry1)[[14]][,5])
#plot(winner_angry1)
#summary(winner_angry1)

# table M.2 model 2
winner_angry2 <- brm(angryprop ~ post*populist_present + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=winners,
                     warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(winner_angry2)[[14]][,5])
#plot(winner_angry2)
#summary(winner_angry2)




#################
# losers + love #
#################

# table M.1 model 3
loser_love1 <- brm(loveprop ~ post
                   + populist_present + polarization
                   + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                   + Xl + Xr 
                   + (1 + Xl + Xr | party_election),
                   data=losers,
                   warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(loser_love1)[[14]][,5])
#plot(loser_love1)
#summary(loser_love1)

# table M.2 model 3
loser_love2 <- brm(loveprop ~ post*populist_present + post*polarization
                   + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                   + Xl + Xr 
                   + (1 + Xl + Xr | party_election),
                   data=losers,
                   warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(loser_love2)[[14]][,5])
#plot(loser_love2)
#summary(loser_love2)




##################
# losers + angry #
##################

# table M.1 model 4
loser_angry1 <- brm(angryprop ~ post
                    + populist_present + polarization
                    + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                    + Xl + Xr 
                    + (1 + Xl + Xr | party_election),
                    data=losers,
                    warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(loser_angry1)[[14]][,5])
#plot(loser_angry1)
#summary(loser_angry1)

# table M.2 model 4
loser_angry2 <- brm(angryprop ~ post*populist_present + post*polarization
                    + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                    + Xl + Xr 
                    + (1 + Xl + Xr | party_election),
                    data=losers,
                    warmup=1000, iter=3000, chains=3, seed=1234)
max(summary(loser_angry2)[[14]][,5])
#plot(loser_angry2)
#summary(loser_angry2)

#save(winner_love1, winner_love2,
#     winner_angry1, winner_angry2,
#     loser_love1, loser_love2,
#     loser_angry1, loser_angry2,
#     file="Stan_15Days.RData")




##################
# result summary #
##################

load("PSRM Replication Files/Stan_15Days.RData")

summary(winner_love1)
summary(winner_angry1)
summary(loser_love1)
summary(loser_angry1)

round(posterior_interval(winner_angry1, "post", 0.9), 2)
round(posterior_interval(loser_angry1, "post", 0.9), 2)

round(mean(winner_love1$data[which(winner_love1$data[,2] == 0),1]), 2)
round(mean(winner_angry1$data[which(winner_love1$data[,2] == 0),1]), 2)
round(mean(loser_angry1$data[which(loser_angry1$data[,2] == 0),1]), 2)

# priors
winner_love1[[3]]
winner_love2[[3]]

winner_angry1[[3]]
winner_angry2[[3]]

loser_love1[[3]]
loser_love2[[3]]

loser_angry1[[3]]
loser_angry2[[3]]


# unconditional effects figure
unc <- data.frame(label=c("Model 1: Winner + Love", "Model 2: Winner + Angry",
                          "Model 3: Loser + Love", "Model 4: Loser + Angry"),
                  est=NA,
                  lwr =NA,
                  upr=NA)
unc[1,2:4] <- summary(winner_love1)[[14]][2,c(1,3,4)]
unc[2,2:4] <- summary(winner_angry1)[[14]][2,c(1,3,4)]
unc[3,2:4] <- summary(loser_love1)[[14]][2,c(1,3,4)]
unc[4,2:4] <- summary(loser_angry1)[[14]][2,c(1,3,4)]
unc$label <- factor(unc$label, levels=rev(unique(unc$label)))

# figure 4
ggplot(data = unc, 
       aes(x = est, y = label)) + 
  geom_errorbar(data = unc, aes(xmin = lwr, xmax = upr),
                lwd=0.8, width = 0.1) +
  geom_vline(xintercept=0, lwd=0.5, lty=2) +
  geom_text(aes(x = est, y = as.numeric(label) + 0.1, label = round(est, 2))) +
  geom_point(cex=4, pch=19) +
  xlab("Posterior Mean of Post Election and 95% Credible Interval") +
  ylab("") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.text.y = element_text(size = 13),
        axis.text=element_text(size=12),
        panel.spacing = unit(1, "lines"))

#ggsave("unconditional.pdf", width = 10, height = 8)


# interaction effects figure
getStdPostCI <- function(model, var1, var2, intr, label){
  coeftab <- summary(model)[[14]][,c(1,3,4)]
  coeftab <- coeftab[intr,]
  
  temp <- data.frame(label=label,
                     est=coeftab[1],
                     lwr=coeftab[2],
                     upr=coeftab[3])
  colnames(temp) <- c("label", "est", "lwr", "upr")
  return(temp)
}

citab <- rbind(getStdPostCI(winner_love2, 2, 3, 12, "Model 1: Winner + Love"),
               getStdPostCI(winner_love2, 2, 4, 13, "Model 1: Winner + Love"),
               getStdPostCI(winner_angry2, 2, 3, 12, "Model 2: Winner + Angry"),
               getStdPostCI(winner_angry2, 2, 4, 13, "Model 2: Winner + Angry"),
               getStdPostCI(loser_love2, 2, 3, 12, "Model 3: Loser + Love"),
               getStdPostCI(loser_love2, 2, 4, 13, "Model 3: Loser + Love"),
               getStdPostCI(loser_angry2, 2, 3, 12, "Model 4: Loser + Angry"),
               getStdPostCI(loser_angry2, 2, 4, 13, "Model 4: Loser + Angry"))

citab$intr <- rep(c("(a) Post Election × Populist Involvement", "(b) Post Election × Polarization"),
                  times=4)

citab$label <- factor(citab$label, levels=rev(unique(citab$label)))
citab$intr <- factor(citab$intr, levels=unique(citab$intr))

# figure 5
ggplot(data = citab, 
       aes(x = est, y = label)) + 
  facet_wrap(. ~ intr, ncol = 2, scales = "free_x") +
  geom_errorbar(data = citab, aes(xmin = lwr, xmax = upr),
                lwd=0.8, width = 0.1) +
  geom_vline(xintercept=0, lwd=0.5, lty=2) +
  geom_text(aes(x = est, y = as.numeric(label) + 0.15, label = round(est, 2))) +
  geom_point(cex=4, pch=19) +
  xlab("Posterior Mean and 95% Credible Interval of Interaction Term") +
  ylab("") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 14),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 13),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        panel.spacing = unit(1, "lines"))

#ggsave("interaction.pdf", width=13, height=6)


# polarization marginal effects
getME <- function(model, variable, label){
  posteriors <- posterior_samples(model, "b_post")
  
  est <- data.frame(t(sapply(1:length(variable), function(x){
    post_est <- unlist(posteriors[,1] + posteriors[,2] + posteriors[,3] * variable[x])
    return(c(variable[x], 
             mean(post_est),
             quantile(post_est, 0.025),
             quantile(post_est, 0.975)))
  })))
  
  colnames(est) <- c("x", "est", "lwr", "upr")
  
  est$label <- label
  
  return(est)
}

range(president_df$polarization, na.rm=TRUE)
polarization_sim <- seq(0, 11, length=100)

metab <- rbind(getME(winner_love2, polarization_sim, "(a) Winner + Love"),
               getME(winner_angry2, polarization_sim, "(b) Winner + Angry"),
               getME(loser_love2, polarization_sim, "(c) Loser + Love"),
               getME(loser_angry2, polarization_sim, "(d) Loser + Angry"))

metab$label <- factor(metab$label, levels=unique(metab$label))

# figure 6
ggplot(data = metab, 
       aes(x = x, y = est)) + 
  geom_ribbon(data = metab, aes(ymin = lwr, ymax = upr), alpha=0.3) +
  geom_hline(yintercept=0, lwd=0.5, lty=2) +
  geom_line(lwd=1.5) +
  facet_wrap(. ~ label, ncol = 4) +
  xlab("Polarization") +
  ylab("Marginal Effect of Post Election") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 14),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16))

#ggsave("marginaleffect.pdf", width=12, height=6)
